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ABSTRACT 


Two new methods for obtaining control and stability deriv- 
atives from observed flight data are developed. The first method 
Is based on a quasilinearization procedure and Is applicable In 
parameter identification problems where the plant Is modeled by 
a system of linear differential equations, and noisy measurements 
of state and control variables are available. Computationally, 
this method Is equivalent to a modification of the Newton-Raphson 
method. The second, a “directed random search* method Is based 
on a concept called evolutionary programming, and Is also applica- 
ble for nonlinear problems. Using X - 1 5 flight test data, the 
two methods are compared and stability and control derivatives 
for the lateral motion of the X - 1 S are given. 
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TWO NEW METHODS FOR OBTAINING STABILITY DERIVATIVES 
FROM FLIGHT TEST DATA 


By George H Burgm 
Decision Science, Inc 

SUMMARY 

Two new methods for calculating stability and control deriv- 
atives from flight test data are developed Digital computer 
programs for both methods have been written and tested out with 
actual flight test data delivered by NASA Flight Research Center, 
Fdwards 

The first method is applicable whenever the system Is repre- 
sented by a set .of linear differential equations and the error 
criterion Is that of minimizing the squared error inteqral It 
uses parameter sensitivity functions (gradients) to obtain an 
algorithm which minimizes the error function This method shows 
very good convergence and can be extended to certain nonlinear 
differential equations It is shown that a modification of the 
Newton-Raphson method will result in the same computational 
a Igori thm 

The second method is a direct search method in the space of 
the unknown system parameters The search proceeds by chanqlng 
one parameter at a time, proqresslnq stepwise to points In the 
parameter space which yield lower and lower error function values 
The most likely successful next step is determined by a finite- 
state machine, which, in Itself, is obtained by a search procedure 
called evolutionary programming 

Obtaining stability derivatives from flight test data is a 
special case of the process parameter identification problem and 
modifications of these two methods are applicable whenever the 
problem of obtaining system parameters from measured (and there- 
fore noisy) state and control variables has to be solved 

The analysis of flight test data of an X-15 flight shows 
that the first method is computationally more efficient than 
the second one Evolutionary Drogramming may have advantages 
in the case of process identification problems dealinq with 
nonlinear differential equations or with nonquadratic error 
functions The first method peimits obtaining error estimates 
for the unknown parameters and the experiments with the qlven 
X-15 flight test data indicate that these error estimates are 
of red 1 1 s t ic size 

It is suggested that these methods be extended to attack 
the problem of obtaining stability derivatives of airplanes with 
non - n eg 1 1 i Labi e nonlinear properties 
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SYMBOLS 

coefficient matrix for the state variables 
(stability matrix) 

element of A 

coefficient matrix for the control variables 
(control matrix) 

element of B 

Inverse of the normal equations 

1 / all dimensional stability derivative 

"17 \ 3 P/ Parameter 

1 I dimensional control derivative 

Tx l 34^1 parameter 


L r ’ L b > L 6R' 


N p> H r • n b* 

n 4 a- V y , 


dimensional stability or control derivative 
parameters analoguos to the preceding defin- 
ition 


F 


linear vector function for the derivatives 
of the state variables 


II 


J 


M 

m 


3z partial derivative of the cost function 
3a j with respect to the j-th unknown par- 

ameter 

(M+l ) = number of observed data points 
number of control variables 


n 


number of state variables 


P 

r 

t 


roll rate, radians/second 
yaw rate, radians/second 
time, seconds 
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SYMBOLS (concluded) 


T 

U 

u 

w 

x 

X 

2 

II 


'‘J 


ft 

l 


5 A 

.1 t 

60 

A 1k 


T 


o 

k 


observation period, seconds 

total number of unknown systems parameters 

control variable 

weiqhtinq factor 

state variable, calculated 

state variable, observed 

cost function 

angle of attack, radians 

j-th unknown system's parameter 

sideslip anqle, radians 

(in connection with evolutionary proqrammlnq 
chanqe in the cost function) 

aileron deflection, radians 

rudder deflection, radians 

dummy control variable 

Kronecker delta 

standard deviation 
bank angle, radians 

Superscripts 

transpose ot vector or matrix 

nominal condition, reference trajectory 

trajectory obtained in k-tli iteration 

a dot denotes the tune derivative 
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INTRODUCTION 


Detei initiation of stability and control derivatives from 
fliqht test data is a special problem from the more qeneral 
area of systems identification, a field which has received 
a great deal of consideration during the last few years A 
short summary and additional references to earlier methods of 
determining stability derivatives from fliqht test data can 
be found ip a recent paper by Lawrence Taylor (ref 1) and in 
an article by Peter Younq (ref 2) 

The work described in this report was performed under 
contract NAS4-1280 bv Decision Science, Inc in San Oieqo 
for the NASA Flight Research Center at Edwards, California 
The purpose of this contract was to use a new computer 
technique, called evolutionary programming (ref 3) for the 
determination of stability derivatives and to compare the 
computational efficiency of this method with the efficiency 
of other, more analytical, methods The analysis of the 
possible mathematical techniques showed that a very efficient 
computational algorithm can be derived by either uslnq a 
quasil inearization technique or a modification of the class- 
ical Newton-Raphson method In adduton to Its efficiency, 
this method nermits estimates of the variances of the 
calculated stability and control derivatives. 

Both methods have been implemented in diqital computer 
proqrams (for a CDC 3600 and an IRI1 360/40 computer) and ob- 
served X - 1 5 flight data have been analysed The assistance 
of Lawrence Taylor and Harriet Smith, both of NASA Fliqht 
Research Center, in conducting this wort is greatly appreciated 


Statement of the Problem 


In the problem of parameter identification of a linear 
system, the process is assumed to be qoverned by the linear 
matrix differential equation 

x - Ax ♦ Ru (1) 

where x is a vector of n state variables and u a 
vector of m control variables, A a coefficient matrix 
(stability derivatives) with dimension (n x n) and B 
a coefficient matrix (control derivatives) with dimension 
( n x m) It is assumed that some or all of the elements of 


4 and B are unknown and to be determined. Unless otherwise 
stated, these elements will be assumed to be time invariant 

In order to obtain estimates of the elements of A and 
B. measured time' histories of the state variables x and 
the control variables u, and possibly of their time deriv- 
atives, are qiven The problem consists now In findinq 
a^'s and b ,j’ s which will, when substituted Into the 

differential equations, match the observed time histories 
if the differential equations are solved with the correct 
initial conditions and the measured control variables u(t) 

In the case of the determination of stability and control 
derivatives from flight data, it can be assumed that the time 
histories have been obtained in carefully planned experiments 
The accuracy and the confidence limits of the calculated para- 
meters depend a great deal on the proper choice of the forcing 
functions, which are, for the lateral motion, rudder and 
aileron deflection These two functions not only have to 
be linearly independent but must have am amplitude which 
compromises between a large value (therefore obtaining a 
good signal to noise ratio in the measured data) and a 
small value (therefore obtaining a motion with negligible 
nonlinear effects) (ref 4) 

If a time history is generated by a process which Is 
governed exactly by a linear differential equation of form 
(1), and if no noise or errors are introduced in the measure- 
ment of x(t). x(t) and u(t), the n(n ♦ m) unknowns can 
be found by formulating n(n ♦ in) linearly independent 
equations using observed values at (n + m) different time 
points, and then solving this'exact system of linear equations 
for the unknown parameters 

It values for the derivatives of the state variables 
are not obtainable by measurement, yet the state variables 
themselves are measured with sufficient accuracy, numerical 
differentiation of the "state variables gives estimates of 
the x(t) and it is, in principle, again possible to 
formulate as many linear equations as there are unknowns. 

However, since the differentiation is a process which 
introduces noise, more reliable results can be obtained by 
setting up more lineai equations than there are unknowns 
and then to solve these equations by a least square procedure. 
A comprehensive summary of least squares methods can be found 
in reference 5, and a short summary of the computational 
procedure lor obtaining least squares estimates of the 
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stability and control derivatives is given in the next section 
j* tins ieport. The least squares technique is used in both 
tne methods described here to obtain a set of starting values 
for the unknown- parameters . 


So far, it has been assumed that time histories obtained 
from processes which can be exactly described by a linear 
system of differential eauatlons of form (1) are analysed 
In most practical situations, however, equation (1) is only a 
first order approximation to a really nonlinear process 
This, of course, is exactly the case if the lateral motion 
of an airplane is approximated by a system of linear differ- 
ential equations Already the existence of purely lateral 
motion is a simplifying assumption, the exact form of the 
equations of motion for airplanes shows that there is 
coupling between the lateral and longitudinal motion, as 
can be seen for fnstance in reference fi. Purely lateral 
motion described by a linear system of differential equations 
also neglects the product terms between roll and yaw and the 
nonlinear aerodynamic forces produced by the control surface 
deflections. 

An attempt to estimate the error bounds on the calcu- 
lated stability and control derivatives must take into 
account the errors introduced by noisy measurement and the 
error in using a linear model of a nonlinear process 


The Linearized Equations of Lateral Motion 

The following system of differential equations is used 
in determining stability and control derivatives for the 
lateral motion (reference 10) 


P = L p p ♦ L r r t L e P 

+ l 4 a'' a + l «r 

+ 

r =t N p P + N f r + N fi i 

+ n 4 a* A * n «r sR 

+ N 

0 «= P - r * Y b 0 + 

V 

+ Y 



( 2 ) 


T’'e last LOluinn of constants are multiplied by a constant 
con**nl force of constant magnitude one These "dummy 
derivatives 41 allow for compensation of drift of the null 
point of the measuring devices 

The A and B matrices follow 


A 


L P 

L r 

L e 

0 

M 

N 

N 

0 

P 

r 

B 


(1 

-i 

Y 




B 

♦ 

1 

0 

0 

0 


B 

l «a 

L aR 

L o 

N 6A 

N 6R 

N o 

0 

0 

Y o 

0 

0 

0 


The element a jl repiesents the angle of attack of the 

airplane during the flight test When performing a lateral 
maneuver, the pilot attempts to keep the angle of attack 
constant Since n can be measured, it is not considered 
as one of the unknours of the elements of the A matrix, 
in the calculation, tie actually (and slightly varying) 
measured value of the angle of attack is used 

The value of Y Q is also measured and therefore a 
known quantity so that there are 14 unknowns to be determined 
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THE LEAST SQUARES METHOD 


The least squares method is based on the fact that the 
differential equations ( 2 ) are valid at every time instant. 
Substituting measured values of p, r, a and 4 and of 
their derivatives into the differential equations permits 
forming simultaneous linear algebraic equations Setting 
up more linear equations than there are unknowns yields 
an overdetermined system of linear equations, the solution 
of which gives the estimates o‘ the unknown parameters 

In the general case where the model of the system is 
given by a differential equation of form ( 1 ). the following 
four steps art required' to obtain least square estimates 

Step 1 In this preliminary step, obtain numerical 
approximations to those derivatives of the state variables 
which are not available from-measurements ' If the observed 
values of the state variables are not too contaminated by 
noise, the following approximation can be sufficient. 

x i^ t k^ 85 [ x i^ t k + i* ‘ x i* t k-i*] / ^ t k + l " l k-l* 

If the observed values are unequally spaced or are very 
neisy, tnqher order approximations may have to be used 
Gcod results can be obtained by smoothing and interpolation 
w-ith spline functions (ref 7) 

Step 2 Number the unknowns in consecutive order, for 
instance 

a 1 1 = “1 


a nn "(n»n) 
b " “(n.n) r 1 

b nm rt (n«n) + (n*m) 


Step 3 Each time instant, at which measured values 
of the 'state variable.- i = l n and the control 

variables u^ i=l m and values of the derivatives of the 
state variables x^ i=l n are available allows (by 

substitution into the differential equations) the formulation 
of n linear equations, in which the are the coefficients 

of the unknown parameters n k , k = l n*n, the Uj are the 
coefficients of the unknowns k=(n*n)+l n*(n+m) The 

values of the derivatives of the state variables form the 
right hand sides of these equations The form of these 
equations is illustrated for a specific case with three 
state variables and one control variable in equation 3 It 
is clear from the form of these equations that they can be 
separated into n independent systems, one separate system 
for every state var-atlc, and theiefore for the unknown 
parameters of one row of the A and B matrix 

Step 4 The overdetermi ned system can now be solved 
by the classical least squares method and the estimates of 
the unknown coefficients as well as their variances can be 
detei mined 

If some of the elements of the A or 8 matrix are 
assumed to be known, the overdetermined system of linear 
equations is obtained by subtracting the products of the 
qiven parameters with their corresponding coefficients from 
the right hand side 

in order to show how the variances are calculated, write 
the overdetermined system of equations in the form 

A n - b (4) 

(one such system for each state variable) 

Assume a vector , for which 

A a 1 = c 

and denote the difference vector between b and c with v 
v - b - c T 

The least squares solution requires v 1 v = minimum and it is 
obtained by solving 

A T A.," = A T b (5) 
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Equations (5) are called the normal equations, and their 
solution can be expressed as 


(A T A) 


A T b 


The variance of the total fit error is given by 

' - yTy 

(M+l )-n 


( 6 ) 


(7) 


where M+l is the total number of equations (see eq (3)) and 
n is the number of unknowns. The quantity M+l-n is called 
the degree of freedom 


The standard deviation for the i-th unknown parameter can 
be calculated as 



( 8 ) 


where c^ is the i-th main diagonal element of the inverse 
of the normal equation 

C - (A T A> - 1 

The value of v T v can be obtained in two ways Obviously, 
by backsubsti tuting the solutions found for the n" into 
equation (4), then summing the squares of the residuals 

O T 

- Cj) yields the scalar quantity v v Computationally 
more efficient is the following way 

v T v = b T b - ( A T b ) u° 

The correctness of tins expression can be seen as follows 

v T v = ( A, i ’ - b ) T ( An" - b ) 

= (,. T A T -b T ) ( A.i " - b) 

- ii 1 ( A 7 A • - A T b) b T An° + b T b 

The first term vanishes due to the normal equations (5) and 
therefore 

v T v b T h - (A T b) T .. (9) 
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> t-^e tne elements A*b are already known as the right 
a' w sides of the normal equations, this calculation 
requires only the calculation of two scalar products 

Finally, a remark about linear dependencies seems 
in order Consider first linear dependencies between 
individual equations (row dependencies) As long as the 
total number of equations minus the number of linear de- 
pendencies is greater than the number of unknowns, the 
linear dependencies are irrelevant for obtaining parameter 
estimates Consider now linear dependencies between state 
or control variables (column dependencies) If the mathe- 
matical model expressed by equation (1) is adequate, linear 
dependencies between the state variables are Impossible 
However, linear dependencies between control variables 
are possible and if two or more control variables are 
linearly dependent over the entire period of observation, 
two or mote columns hi the overdetermined system, and 
therefore also in the normal equations will be linearly 
dependent In this case, it is not possible to determine 
all the unknown parameters in the B matrix, only ratios 
between parameters can be calculated 

In a process loenti *'<‘ation problem, where the time 
histories are obtained by performing carefully planned 
experiments, linear dependencies between control variables 
can always be avoided For the determination of the stabil- 
ity and control derivatives of the lateral motion of an 
airplane, linearly independent aileron and rudder deflections 
will quarantee a nonsingular coefficient matrix of the 
normal equations — 
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THC METHOD OF QUAS I L! NEAR 1 ZAT I ON 


The method developed in this section results in an 
algorithm which is computationally equal to the one described 
by Taylor in his paper “A Modified Newton-Raphson Method for 
Oetermimng Stability Derivatives From Flight Data" (ref 1) 

It is interesting to note that the same procedure for cal- 
culating stability derivatives can be obtained by two quite 
different approaches 

The basic idea is to find coefficients in the A and 
B matrix which permit fitting the observed time histories 
(which are assumed to be solutions to the differential 
equation (1)) As criterion of fit the integral of the 
weighted squared differences between calculated and observed 
time histories is chosen The following cost function is 
therefore defined 

T 

2 3 £ w, /fx,(t) - x,U)] dt (10) 

, = 1 0 

I* reliable measurements of fhe derivatives of the state 
vaiiables are available, a different possible cost function 
could be taken as 

T 

z * £ **, / [ *,(»)-* rt)]’’ dt 
1 = 1 o L 

f 

* L «, J [ *,<t) - *i(t>]' dt on 

1 = ] 0 

It is also possible to include only certain derivatives 
in the second sum. eg p and r It seems that the choice 
of the tost function deserves additional attention 

As weighting factors the inverse of the root mean square 
of the observed state vaiiable qives a reasonable balance 
between the four observed state variables Different choices 
for the weighting factnr r aie. or course, possible and may 
take the telative accuracies of the measurements Into consider- 
ation (lor instance, roll and yaw rates can often be measured 
more accurately than bank and sideslip angle ) 


13 


In the following derivation, a cost function of form 
. 'O' is assumed and for simplicity of notation, the weight 
ravtois ate all assumed to be one 

As a preliminary step in deriving the algorithm It 
is shown how the sensitivity functions (here the partial 
derivatives of the state variables with respect to the 
unknown parameters) can be calculated Rewrite equation 
(1) in the following form 

x = F(x. u, o, t) (12) 

where f is a vector function and all unknown parameters 
are combined Into a row vector a, such that 



Differentiate equation (12) with respect to some a. 
say 1 ,^ , and write the 1-th component of the vectu ( n 

assulned to be Independent of time) 


3 F f (x, u, o, t) 
' ,a J 


y ♦ e ^ 

3x k 3a j k = i * u k 3a j 


+ 


(13) 


Since the control variables arc independent of a, the second 
sum vanishes Under the usual assumption that the second 
paitial derivatives are continuous, we can interchange the 
order of differentiation and obtain 


X 1 

n 

E 

,F , 

DX 

3 F 

) 

T^i 

lx k 

.1 u 

3n 


(14) 
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Equation (14) is a linear differential equation for the 
influence coefficient 

DX, 

JOj 

Specifically, for a system of form (1) with four state 
variables and three control variables, the sensitivity 
equations can be written as 





1 = 1 .4 

k-1 , . 4 
J-1 .3 


(15) 


where is the Kronecker delta 

There are 108 linear differential equations, which 
can be solved simultaneously with the four equations (1). 
Therefore, a system of 112 differential equations Is 
obtained The initial conditions of the Influence coef- 
f iclents 



are all zero because the parameters are Independent of 
the initial conditions of the state and control variables 
Equation (lb) indicates how the sensitivity functions are 
obtained as solutions of a linem system of differential 
equa t Ions 

It is now shown how these influence functions can be 
used to obtain collections to the unknown coefficients 
Obtain a first appi ox .mat ion to the unknown (constant) 
coefficients bv applyinq the least squares technique as 
described in the preceding section Then solve the differ 
cntial equations for the state variables together with 
the equations for thp para m ter influence coefficients 
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At each point at which observed values of the state 
variables are available, expand around the reference' 

-~oidt obtained with the present values of the coefficients 
the solution of the state variable in a Taylor series 
as a function of the unknown parameter corrections, 1 e.. 


*.,(» + An, t) = x 


n "1" + 

, U.o - £ 


( u , t ) An 


* higher order terras 

In the above expression, consider x (a+Ao,t) as the desired 

, o 

value (equal to the observed value of x^ft), x^u.t) the 

value of the presently computed reference solution and the 
summation as the desired correction This clearly gives a 
linear equation for the A<,j . For each state variable, and 

for each point t, one such equation is obtained When 
the reference solution is carried out over the entire ob- 
servation interval, n(M+l) linear equations can be formulated 
and solved by a least square method This solution yields 
corrections to the present values of the parameters Adding 
these corrections to the parameters will give a new reference 
trajectory, closer to the one of the observed data This 
process can be repeated until the corrections become neglig- 
ible. 



One way of lookinq at the problem of obtaining corrections 
requires that the n integrands in 




( a , t ) Aa 


(U is the total numbei of unknowns) 
vanish In other words, try to satisfy the following equation 

H ( a , t ) An = 

j = l 3 *'j J 

Formulating these equations for t = t 0> yields an 

ovcrdeterunned system for which the normal equations are of 
the form 


x - x °(o,t) 
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and an IBf .160/40 computer Four to five Iterations were 
lequued to get the corrections to about 1/10,000 of the 
\alue of the coefficient Using approximately 120 time 
points required approximately one minute computer time 
on the CDC 3600 to calculate 15 unknown parameters with 
5 iterations The differential equations are solved by 
a fourth order Runge Kutta method 

Estimates of the variance of the unknown parameters 
are obtained in the following way 

The' variance of the total fit error can be expressed as 
O' - £ ", •£[*,<»„> - ?,<»„)] '/(M + l-U) 

TT, v=o* 

and estimates of the variances of the individual parameters 
are calculated as 


(19) 


( 20 ) 


where Cjj is the j-th main diagonal element of the Inverse 
of the normal equations (18) 
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V * 


? X 




1 (t )*.. ♦ £ £^(t ) £-<t ) A 


, 5TZ-. " a« ' v' 3a ' v' 

i - 1 v -0 ? \ 


a + 


■ E E [ s(‘ w )-« i (t„)] 

1=) v-o J 1 


'I M 


i M 


fla ♦ 
2 


E E ’u ) -« ( t ia. ♦ E EPit ) 3 

i- i tt v 3a v i fz* TTT. Da v 

1=1 V P 2 1 1 = 1 v = 0 2 

■ 5, s„f !r ' I, . : £"•> 




" M 


E E^(t„)^(t j*. ♦ 

Vi „. 0 U u 3a - v ? 


(18) 


The solution to these normal eauations yields the corrections 
\, tj , which are added to the old values of the o j and then 

a solution of the differential equations for the state varia- 
bles and sensitivity functions using these new parameters 
is performed This process is iterated until the corrections 
become negligible 

A computer program which allows up to 15 parameters 
assumed to be unknown was written and run both on a CDC 3600 
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NEWTON-RAPHSON'S METHOD AND ITS MODIFICATION 


This section shows how the same computational algorithm 
which is given in equation (18) can be obtained by a modi- 
fication of the Newton-Raphson method and that It might be 
worthwhile to program the unmodified Newton-Raphson method 
After a short general exnositlon of how the Newton-Raphson 
method can be used In optimization problems. It Is shown 
how the necessary partial second derivatives of the state 
variables with respect to the unknown parameters can be 
calculated In a manner similar to the one applied for ob- 
taining the first order sensitivity coefficients It Is 
then shown how the second order partial derivatives of 
the cost, function with respect to the unknown parameters 
are obtained a'nd used to minimize the cost function 

First, consider the problem of minimizing a function 
of several Independent variables with no constraints 
Let 


Z = F(<*^ a n ) ■ F(n) 


where « denotes the vector with elements a , a„ 

l u 

A necessary condition that z has a local minimum Is 


, = 0 for all j 

"“J 

Since Newton's method is really a procedure to find zeros 
(not extrema) of functions, it is used to find values of 
« which will satisfy the above condition Let 


V‘‘> 


3Z 

'"j 


and expand 11^1*) around some point a* into a Taylor 
series 


HjU 


■'‘.I 


V : ') 


u 

£ 


3H, 

r„J (a ° )A “k 


+ higher order terms 


(81) 


(? 2 ) 
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Meolecting the higher order terms and requiring that 
(a° + Au ) a 0 

yields a system of U linear equations for the U unknowns 



( «" )Ao k 


v»-> 


J=1 u (23) 


Newton's method consists in solvinq the above syster of 
linear equations. Since the higher order terms have been 
neqlected, H J (u n + An) will not be exactly zero but, if the 

starting point was close enough to the zero of Hj , 

|H ( a ♦ iu)| will be smaller than |h (a”)| The process 

ts repeated and converges to the zero of with quadratic 
convergence Remembering now that is the first partial 

derivative of the function to be minimized with respect to 
the j-th unknown parameter, the elements in the coefficient 
matrix in the above equation are the second partial deriv- 
atives of the function to be minimized with respect to the 
unknown parameters In order to obtain the second order 
partial derivatives of the cost function, consider first 
the procedure to obtain second partial derivatives of the 
state variables Write the differential equation which 
governs the state variables in the following form 



The (vector) function G may be nonlinear Then 



(24) 


20 


Assuming the control variables independent of a, the 
second term vanishes and since the a are assumed to 
be time independent, the third term drops out Inter 
changing the order of differentiation, we obtain the 
well known result 


d / ax i \ y aG i ”_k + 3 _ G i 

dt \n.ij / ax^ acij aa^ 


Differentiating again and immediately changing the order 
of differentiation on the left-hand side gives 



a 2 6 i 

3x u 

3x k 

3 j G, 

4 * 

3x k 

T5T 3xT 
u * 

30 t 

3 °J 

3 V x k 

5a 

J 


+ a2 x k i’S + 3?G i \ + 3?6 1 (25) 

an^a^ ^ x |(, 3 a j 3x k 3o t / aa j 3a t 

This system of differential equations, together with the 
differential equations for the state variables and the first 
order sensitivity functions allows the calculation of 
the second order partial derivatives It may be emphasized 
again that the above deri vat ionmade no assumption about 
lineaiity of the functions G^fu, x, u. t). 

The total number of second partial derivatives and 
therefore of differential equations of form (25) required 

in the Newton-Raphson method is nil 2 

for the special case where the functions G 1 are 
linear and of the form 


a U X J <t) 


■evcral simplifications can be made 


/ 
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= 0 


for all s and t 



The differential equations for the second partial 
dnivdtives aie 


■ ll (•'.l’*.,) 

£ *■* 


k i .i, v ♦ h 

ld st' >a uv iu ’ d st ’ 


rt,k 




b \x 3tr yz 


(26) 


In all of the above expressions, i, s, t, u, v, w, and y 
run from 1 to n, x, and z from 1 to m 

It is appropriate to mate a remark about the order of 
magnitude of the task of calculating the second order partial 
derivatives needed In the Newton-Raphson method Assume 
n a 4 and m * 3 and assume 7 unknown parameters in the 
A and 7 unknowns in the B matrix Due to the symmetry 


i)o j 3cijj 


a total of 4*14*15/2 • 420 

second partial derivatives are required, which means that 
a system consisting of 480 (linear) differential equations 
(4 state variables, 56 first order partial derivatives and 
420 second order partial derivatives) has to be solved Con 
sidering the simple form of the riqht hand sides of these 
equations, it is feasible to solve them on a digital computer 
with typically 32,000 words of core storage 

The above derivation showed how the second partial 
deiivatives of the state variables with respect to the un- 
known parameters can be found The necessary groundwork 
is now laid to consider the problem of determining stability 
derivatives using Newton's method 
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Let 


2(0 

■ 2 Z /fxjt.n) - x ( t ) 1 dt 

IS| ft ■* 

(27) 

Identify 

1 1 0 

n 7 


iz 

'j 

2 iff \ x , (t - ; > - *, 10 ] dt 

0 *i 

(28) 


with the function H j(") of the first paragraph 

Let again <■' be a point in the parameter space close to 
a local minimum of Then the n linear equations for 

obtaining the are (see equation 23) 

U (i T 

5a] « = *•*) «}^« k 

n T 

= - ' JJ /f*i<S".t) - ?,(t)]^’ (S".t) dt 

° J (29) 

J = 1 U 

Cairyinq out the <1 i f f ei en 1 1 a 1 1 on under the sum and inteqral 
sign and dropping the constant factor ? gives 

U " V« 

k? ( '‘’ t) + [ X i ( * A) - 

I! T t 

dt 'k ' £ /f x i ( " ’ l) -*■, < *■> ] (« .t) dt 

1 1 I) 1 1 J 


J = 1 U 


Z4 


\**c*ing now 



(o.t) , 


the coefficients in the matrix for the linear equations can 
be obtained by an integration of known functions Also the 
elements of the right hand side vector can be obtained by. an 
mtegi at ion of known functions 

The Modified Newton-Raphson Method 

If the integrals in equation (30) are approximated by 
sums and if the second partial derivative term is neglected, 
the following equations are obtained 



(31) 


for j=l U 

Comparison of equation (31) with equations (18) shows 
that the two systems of equations for obtaining the corrections 
I >* . are identical 

This demonstrates that the application of the Newton- 
Raphson method to the minimization of the cost function (10) 

gives the same result as the method of quasi 1 linearization H j 

the second partial derivatives are neglected The possibility j 

of modifying the Newton-Raphson method by neglecting the second I 

partial derivatives was mentioned the first time by Balakrish- 4 

nan in reference 8 
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DIRECT SEARCH BY EVOLUTIONARY PROGRAMMING 


Basic Concepts of Evolutionary Programming 

Evolutionary programming is described in detail in 
the book "Artificial Intelligence Through Simulated Evolution" 
by Fagel , Owens and Walsh (ref 3) In order that a reader 
unfamiliar with the concept of evolutionary programming might 
understand the direct search method described in the next 
section, a short summary of evolutionary programming is 
presented 

Consider first a Moore-machine , a triplet (I.i.f), where 
I represents an input alphabet with a finite number of ele- 
ments, z a set of states and f a transfer function from 
I x ) to the set of states When an element of the input 
alphabet is received by a Moore-machine, it will transfer 
from one state to the next state in accordance with the rule 
laid down by the transfer function f The program used for 
the determination of stability derivatives is capable of 
handling Moore -machines with up to five states and with an 
input alphabet size of 60 Once the input alphabet and the 
number of states are specified, the transfer function f 
is given in tabular form specifying for each state and for 
each input symbol the next state reference 

Evolutionary programming works essentially with finite 
state machines which can be described by a quintuplet 
(l.t.f.O.g), where again I denotes an input alphabet, I 
a set of states, f the next state transfer function, 0 
an output alphabet which may or may not be identical with 
I (in the application for stability derivatives, 0 con- 
tains only 6 different elements), and g an output 
function For a given output alphabet 0 the specification 
of g will uniquely determine a finite-state machine 

Evolutionary programming is a method to find finite- 
state machines which will produce, for a given sequence of 
ir"'ii symbols a sequence of output symbols which will minimize 
a certain cost function Evolutionary programming consists 
essentially of three basic procedures, an environmental 
comparison procedure, a mutation and selection procedure, 
and an output determination procedure 



EXAMPLE OF A FINITE-STATE MACHINE 

3 internal states, 3 symbol input alphabet, 2 symbol output 
alphabet 


Figure 1 shows a finite-state machine with three states, 
three input symbols and two output symbols Assume that the 
finite-state machine is initially in state A and receives 
the input sequence 0020 The following sequence of 

events will then take place 


Present State 
Input Symbol 
Next State 
Output Symbol 


ABB 
0 0 2 

BBC 
0 1 0 


C B C A 
0 2 2 1 
B C A A 
0 0 10 


The output determination algorithm drives the Moore-machine 
with the given sequence of input symbols and determines those 
outputs (to each branch) that will minimize a given cost 
function The outputs are obtained in a deterministic manner 
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T h «* imit.it ion and selection procedure generates an off- 
••.•iiKi In i .iiul mill y pet f orm i ng one of the following mutations 
of the L’lesenl finite-state machine 

Adding a state 

Deleting a state 

Changing a next state reference 

Changing the start state 

When a state is added, some next state references to 
previously existing states are changed randomly to connect 
the new state with the rest of the machine. The newly 
added state obtains as many next state references as there 
are input symbols These next state references are added 
randomly When a state is deleted, all the next state 
references which referred to that state are deleted and 
randomly connect to some other states 

The two remaining mutations are self-explanatory The 
output determination routine determines now the proper output 
symbols to this mutated machine 

A comparison is made whether the parent or the offspring 
obtains a lower value of the cost function and that finite- 
state machine yielding the lower value is kept as a new 
patent ma ch i ne 

As more and more information becomes available (larger 
and larger recall) over which the finite-state machines can 
be exercised, finite-state machines are generated which re- 
flect, with increasing fidelity, the logic of the underlying 
process 


Evolutionary Programming Applied to Function Minimization 

The basic idea of evolutionary programming is to find 
finite-state machines which reflect in some sense the logic 
in the behavior of a system This may be an independent 
system or a system which interacts with its model and whose 
behavior, therefore, is dependent on the evolutionary program 
As an example of the first class, consider the problem of 
finding finite-state machines describing the logic of the 
changes in ocean temperature Clearly, the ocean temperature 
is independent of the logic found by the evolutionary program 
The problem of function optimization is an example of the 
second kind, such a process can be viewed as being a game 
between the evolutionary program and the cost function. The 
behavior of the value of the cost function in the past is now 
dependent on what "moves" the evolutionary program made, this 
is the interaction between the two "players" There is of 
course a clear distinction between the cost function and the 
"values of the cost function” The cost function itself is 
certainly independent of the optimization method used, but 
the "values of the cost function" depend on the path taken 
by the optimization procedure It may be mentioned here that 
Wilde, in his book on optimum seeking methods (ref 9) also talks 
about the "Opening Gambit" and the "End Game”. The purpose 
of the evolutionary program in an optimum seeking procedure 
can be summarized as being the device which indicates which 
sequence of changes in the free parameters will be the most 
promising to reduce the value of the cost function This is 
done presently in the following way 

Each free parameter can be changed in four different 
ways, in a positive or negative direction with either a 
large or small step size Consider now an input alphabet 
( 0 ^ ) consisting of four times as many symbols as there are 

unknown free parameters Each of these symbols represents a 
unique change in one of the parameters An evaluation of the 
cost function using this changed parameter will yield a new 
value of the cost function Designate the change in the cost 

function by , therefore, « z 1 - z y -j Negative values 

of b correspond to improvements in the set of parameters, while 
a positive value means a degradation in the set of parameters 
Clearly, it is desirable to make B as small as possible 
(this corresponds to a large improvement) Define the follow- 


28 


29 


mg intervals 

a 1 < a 2 < o < a 3 < a* 

and define a corresponding output 6 where 


s 

6 

H 

If 

I' 

n 


1 If 

8 < a 

l 


V If 

a i < B 

< 

a? 

1 if 

a? < b 

< 

0 

If 

0 < B 

< 

® 3 

•. if 

a j < it 

< 

at, 

( 1 f 

a„ < 8 




for the given set of parameters there corresponds exactly 
one output symbol p to one input symbol a 


Assume now that there exists a finite past history of 
pairs of { .1 k , h^) Clearly, the game of minimizing the cost 

function with respect to the given parameters consists now 
in choosing an u k + j which will produce a 6 ^ + -| as small 

as possible It Is here the evolutionary program comes into 
play Suppose that there is a finite-state machine which 
will "fit" the sequence (o n> B n ), n = 1, 2 k fit here 

means that if the finite-state machine is driven by the se- 
quence of the a n , it will produce the output sequence 8 n 

Then at the k-th move, that finite-state machine Is In a 
certain state, say All possible inputs B will have a 

unique output associated with them It seems logical to assume 
that, since the finite-state machine was a perfect fit over 
the past, this finite-state machine contains information 
about the outcome associated with any given next input symbol 
Scanning all possible outcomes and searching for the lowest 
possible, the input symbol associated with the lowest output 
symbol can be determined Call the lowest possible output 

symbol 1 k + 1*^ The s y n| hol associated with R P ^ P< * 15 now con- 
sidered to be the evolutionary program's next move Note that 
several different input symbols may produce th-c '-me lowest 
value for the output symbol If this is the case, one among 
all the candidates for producing the lowest output symbol Is 
chosen randomly foi actual usage The parameter change assoc- 
iated with this symbol ^ + j is performed and a new evaluation 

of the payoff function occurs At this point, the two newly 

generated symbols „ k + | and (which corresponds to the 


actual change in the cost function), are added to the list 
of Input/output pairs If the actual output, bJ| p *. was a 3 

or smaller, the game proceeds In the same way as described 
above, the recall being now one event longer If the output 
was 4 or more, the evolutionary program Is considered as 
having made an action error. This case Is described later 
in this section. 


It Is shown above how, in principle, a finite-state 
m-c-hine with a perfect fit over the past history is used 
as an “acting player" It is appropriate at this point to 
add some additional important details First, for any 
reasonable length of the sequence of the past moves, called 
the recall, it Is unlikely that a finite-state machine with 
a perfect fit will be found For this, and other reasons, 
more than one finite-state machine are carried as possible 
players In the evolutionary program, specifically, for the 
problem solved here there are three finite state machines 
which are restricted In size. Machine 1 Is a simple one-state 
machine while machines 2 and 3 can have any number of states 
between 2 and 5. Before a move Is made, the “best" one 
of the three possible machines Is selected as "player". Best 
means that machine with the lowest fit score The fit score 
Is obtained In the following way. Given k pairs of Input/ 
output symbols, B)> a 2 , fl 2 , a k . 6 k , drive the finite- 

state machine with the k inputs and for each move form the 
difference I B m - B a l, where B m is the output predicted by 

the machine and h. the output that actually occured Divide 
k « 

S ] J 

il B n B al by k and define this quantity 


as the fit score If the actual output Is qreater 


than or equal to three (a degradation in the cost function) 
the evolutionary program Is considered to have made an error 
or a 'bad action" has occurred Machines 2 and 3 are now mutated 
Mutation means that with probability 


p l 

p 2 

p 3 

p 4 


a start state Is changed 
a next state reference Is changed 
a state is added 
a state Is deleted 
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where, of course 



If the mutant has aliuady the maximum number of states, 

1 ^ automat Iia 1 I y is set to zero and analogously equals 

zero i< the machine has only 2 states Machine 1 is not 
mutated, because the only possible mutation would be to 
add a state, but the intention is to keep machine 1 a one 
state machine 

The offspring machine is driven over the recall and its 
fit score is evaluated If the fit score Is worse than that 
of the parent machine, the offspring is replaced by the 
parent and another mutation is tried A parameter in the 
program limits the number of trials for obtaining a better 
machine After machines 2 and 3 have been mutated (or at 
least an attempt has been made to mutate them), the machine 
with the best fit score of the three machines is chosen as 
the actor and the game proceeds in its normal way 

Two important details of the evolutionary program have 
not yet been discussed, the setting of the outputs and the 
treatment of the unexercised sta te- trans l t l ons First note 

that the mutations affect only the structure of the finite- 
state machines but not their outputs It is clear that for 
a finite-state machine with a given structure and a given 
start state, there exists at least one setting of the outputs 
which will minimize the fit score Assume first that no 
state tiansition is exercised more than once during the 
recall Then, each one of the exercised s ta te- trans i s t i ons 
is assigned a unique output symbol (namely the one corres- 
ponding to the actual occurred output symbol during the 
transition) and the fit score will obviously be zero, because 

for all j, j = 1 k, n J = More important is the case 

d rn 

where some st a te- trans 1 1 i ons are exercised more than once 
during the recall It is possible that, in order to fit the 
actual data, different symbols would be required each time 
the transition is exercised In such a case, the output 
symbol is sot to some weighted average (rounded to the near- 
est integei) of the desired output symbols Also a new fit 
score is detined, which is equal to the fit score as described 
above divided by the number of times that non-unique state- 
tiansitions have occurred, this is called the normalized fit 
Scot e and it is this normalized fit score on which the choice 


of the acting player is based. The third possibility is 
that a state-transition is never exercised during the 
recall and therefore the output symbol associated with 
this state-transition has no influence on the fit score 
It proved to be advantageous from a programming point of 
view to assign a special symbol to those outputs An 
output of zero now designates an unset output value 

A last point to discuss is the procedure used if 
analyzing the acting player, a move that will produce an 
improvement in the cost function cannot be found, or 
expressed in the alphabet of the evolutionary program, 
if in a given state no input symbol will produce a pre- 
dicted of 3 or better. If this occurs, the evolution- 
ary program is not used to generate a symbol for the next 
move The next move is obtained by scanning the past moves 
and finding the most recent move which produced a B of 
3 or less If at any move, an actual output of 4 or 
greater is generated, the evolutionary program 1$ said 
to have made an action error If this occurs, the next 
move or next moves will not be determined by the evolution- 
ary program, but rather by a subprogram which essentially 
tries out whether a step in the reverse direction is better 
and it keeps trying until it again finds a successful move, 
always restoring the coefficients to their old value after 
an unsuccessful move Details about this subroutine can be 
obtained from the flow chart in Figure 2 This subprogram 
also guarantees that at the end, a local minimum of the 
cost function within the specified levels of changes in the 
parameters has been found, because only after an exhaustive 
unsuccessful search over all possible single changes In 
the parameters is the run terminated A second way to 
terminate the run is by limiting the number of moves After 
termination, the final value of the coefficients are printed 
together with the complete sequence of input/output pairs 
Since the evolutionary program requires some environment in 
the past, initially, to start the program, a separate sub- 
routine generates a prescribed number of input symbols and 
the corresponding output symbols are calculated. This is 
called the initial environment The moves which generate the 
initial environment are performed in a stochastic manner 

Summarized, the features of the present version of the 
evolutionary program for the determination of stability and 
control derivatives are as follows 


(text continued on page 37) 
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FIGURE 2 

HOUCIIART OF I VOLIIT I OHARY TROCRAN FOR FUNCTION MINIMIZATION 
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FIGURE 2 (CONTINUED) 
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differential Equation x = Ax + Bu where 

A Is of dimension (4,4) 
B is of dimension (4,3) 


Cost Function - j 


Z = ]£ w 
1 = 1 

-7 

(x 1 

- 7 i 

)* dt 


Free Parameters a)) 

a ) 7 

a 1 3 

a 2 1 

a Z Z 4 2 3 

4 31 * 33 

b n 

biz 

b t 3 

b 2 i 

b Z 2 b 2 3 

b 3 3 


Number of Finite- 

State Machines - 1 one-state machine 

2 2 to 5 state machines 


'Input Symbols 60 

Output Symbols 6 


Input to the Program 

For each machine Initial configuration 
Initial start state 
Maximum recall length 


Probability distribution for 
the 4 types of mutation 


p, = probability of-changlng 
start state 

p- = probability of changing 
next state reference 

p, - probability of adding a 
state 

p. * probability of deleting a 
state 


Maximum number of moves 
Print Interval . , 

Number of errors allowed before a mutation occurs 
Maximum number of machines tried at a mutation 
For all 60 Input symbols the change in the coefficient assoc- 
iated with this Input symbol 

The interval limits for the determination of the output symbol, 
a, , 1 = 1 . 4. 

The initial values of the 28 elements of the A and 8 matrix. 
The observed time histories 
The integration step size 
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Resul ts 


some key results are listed below for a typical run of 
t hr evolutionary program on the CDC 3600 computer 


Running time 



6 

minutes 

Initial environment length 

50 



Total number 

of 

moves of the evolutionary 




program 



400 



Total number 

of 

function evaluations 

750 



Value of z 


initially 


1 

81 1 



after mit 50 moves 


1 

038 



after 100 moves of Evol 

Pr 

0 

564 



after 200 moves 


0 

491 



after 400 moves 


0 

463 


Input Symbols The 60 symbols of the input alphabet represent 
changes of 


of 

+21, 

the coefficients 

CSJ 

+ 0 2t , 

-0 2% 

a J 

i ' a i / • a i i • 

a . 1 • 

a ,,, 

a r 3 * a 

31 * 3 3 3 


1 • b 1 , • b l 3' 

b > i * 

b„. 

b ?3 . b 

3 3 

1 n 

this order 






Intervals for output symbols 

Output Symbol Corresponds to 


1 


< 

-5 

0 

10- 3 

2 

P 

< 

-5 

0 

10- 4 

3 

H 

< 

0 



4 

e 

< 

1 

0 

10‘ 3 

5 

V 

< 

1 

0 

10~ 2 

6 

p 

> 

1 

0 

eg 

■ 

O 


The initial coefficients were those found by the least squares 
piocedure and are given in the following two matrices 


A 


-0 099*) 

0 595 

-22 48 

0 

0 00641 

0,0668 

1 036 

0 

0 1 147 

-t 

0 02 38 

0 00698 

1 

0 

0 

0 
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B 


12 99 

15 11 

0.361 

0 487 

-1.764 

0.00731 

0 

0 

-0.00272 

0 

0 

0 


Note that In this experiment a 31 was considered as one 
of the unknown parameters. 

The coefficients at the end of the run were 


A 

-0 153 

0.163 

-21.67 

0 

0.00615 

0.0642 

1 108 

0 

1.216 

-1. 

-0.0190 

0.00698 

1 

0 

0 

0 

B 

13 25 

14 53 

0 378 


0.551 

2 06 

-0.00221 


0 

0 

-0.00261 


0 

0 

0 



Note that the first eleven moves of the evolutionary program 
(move 51 through 61) produced all outputs of 1, which Is 
quite remarkable, considering that the longest string of 
"1" in the first fifty moves was only of length 3 (see Table I) 

Although the optimization method using the evolutionary 
program works satisfactorily In its present form, there exist 
possibilities to improve Its performance. A first Improvement 
consists in preventing the evolutionary program from qettlnq 
"trapped" In a long string of input symbols which all produce 
an output 3 (a very slight improvement In the cost function) 
If unlimited computer time were available, these long strings 
of outputs of 3 would be all right, but In the Interest of 
of saving computer time, an attempt should be made to find 
parameter changes which will improve the cost function more 
rapidly. Such values may be found more quickly If after a 
string of outputs 3 with some given fixed length, a random 


(text continued on page 42) 
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*earch procedure similar to the one used to generate the 
initial starting, environment, is used foi j given number 
of moves 

A second improvement would, as the search procedure 
approaches the minimum of the cost function, automatically 
change the magnitude of the changes in the coefficients 
(say reduce them by a factor of 10) and also reduce the 
values of |a ( | throuqh | a ,, | This latter change would 

increase the sensitivity of the evolutionary program to 
changes in the cost function 


EXPERIMENTAL RESULTS AND COMPARISON OF THE TWO METHODS 


Two computer programs were developed, one Implementing 
the quasi 1 Inearization method and the other one the direct 
search method using evolutionary programming 

The first experimental runs with these programs were 
made with data which did not originate from actual flight 
tests, but which were obtained from solving four simultaneous 
differential equations of form (1) on a hybrid computer and 
using the measured and digitized data from these runs. Clear- 
ly, since these data originated from a process described 
exactly by a differential equation of the form considered 
hece, and since the only errors were roundoff errors In the 
digitized data, the coefficients were found quite accurately 
and the observed time histories were matched by the calculated 
time histories with the same accuracy as the originally g-tven 
time history data (three to four significant digits). 

The next case analysed the flight test data of an X-15 
flight. Measured data of p. r, b, a, and of 0 amp -r"Were 
available at 0.025 second Intervals for a total observation 
time of 6 seconds. For the calculation, every second point 
of these time histories was used. 

A first approximation to the unknown coefficients was 
obtained using the least squares method In the experiments 
with the program using the quasi 1 1 near 1 zatl on method, for 
the element a 3 i the observed angle of attack has been used. 

By the least squares method the following parameters and 
estimates of their variances were obtained. 


A 


-0.101*0 on 

0 539*0.213 

-22.43-0 15 

0 

0 0064*0 001 

0 0619*0.019 

1 .036*0.01 

0 

0 114*0 0019 

-1 

-0.058*0.021 

0 00698 

1 

0 

0 

0 


B 


12 99*0 40 

15 15*0.38 

0.359*0.006 

0 498*0 035 

-1 760*0.034 

-0.0074*0.00058 

0 

0 

0.0148*0.0013 

0 

0 

0 
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After five iterations with the quas 1 1 i nearl zat ion 
method, the following coefficients and error estimates 
weie obtained (see also Figure 8 for comparison) 


A 


-0 1 9 1 fO 04 

2 853*0 75 

-24 08*0 24 

0 

0 0041*0 0028 

-0 126*0 06 

0 974*0 025 

0 

u(t) 

-1 

- 0 020t0 035 

0 00698 

1 

0 

0 

0 


B 


14 21+1 45 

19.37+1 72 

0 406+0 025 

0 709+0 128 

-1.951+0 159 

-0 002*0 002 

0 

_ 0 

-0 0012*0 0008 

0 

0 

0 


It was beyond the scope of the work performed under this 
contract to develop methods for obtaining error bounds on 
the calculated stability and control derivatives Nevertheless, 
the methods used to find numerical values for the variances in 
the calculated derivatives seem quite reasonable and they allow 
at least an estimate of the expected relative accuracy of the 
parameters. For instance, looking at the value of a )2 (L p ) 

and its estimated variance in the pure least squares solution 
Indicates that this parameter was determined with very little 
accuracy Indeed, the a | 2 found after the fifth iteration 

differs from the a 12 from the least squares solution by a 

factor of about five, and again, the estimate of the error 
after the fifth iteration is still fairly large On the other 
hand, looking at a 73 (N e ) the relative small variance in the 

least squares solution is an indication that this parameter can 
be determined relatively accurately and the final value of 
a 2l after five iterations differs only about 6? frem the value 
found by the pure least square procedure 

Tiqures 3 through 6 show the observed roll and yaw rates 
and the sideslip and bank angles On the same graphs are shown 
the time histories obtained using the coefficients of the pure 
least squares solution and the trajectories obtained with the 
coefficients after five iterations of the combined qradient- 
least squares method Figure 7 shows the corresponding aileron 
and rudder deflections 
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FIGURE 3. 

OBSERVED AND CALCULATED ROLL RATES OF X-15 FLIGHT TEST. 


a - Observed roll rate 

h - Calculated roll rate using coefficients found with 
pure least squares procedure 
t - Calculated yaw rate usinq coefficients found with 
5 iterations of the quas 1 1 1 nearl zati on method. 


( text continued on page 51) 
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Time, seconds 
FIGURE 4 

OUSIRVID AND CALCULATED YAVi RATES OF X - 1 5 FLIGHT TEST 
a - Observed yaw rate 

l> - Calculated yaw rate using coefficients found with 
pure least squares procedure 
c - Calculated yaw rate using coefficients found with 
S iterations of the quasi t i neari zation method 


Side 



lime, seconds 


FIGURE S. 

OBSERVED AND CALCULATED SIDESLIP ANGLE OF X-15 FLIGHT TEST. 


a - Observed sideslip angle. 

b - Calculated sideslip angle using coefficients found with 
pure least squares procedure. 

c - Calculated sideslip angle using coefficients found with 
5 iterations of thequasillnearizatlon method. 







radians 




FIGURE 8 
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The experimental results of the evolutionary program 
were discussed in the section "Oirect Search By Evolutionary 
Programming". 

Figure 8 compares the values of the coefficients obtained 
by the least squares method with those obtained by quasillnear- 
ization, as well as the corresponding estimates of the var- 
iances The fact that the estimated variance in the least 
squares method is smaller than in the quas i 1 1 near 1 zatlon method 
should not be taken as an indication that the least squares 
coefficients are closer to the true values than the ones 
obtained by the quasi 1 inearization method It seems that the 
parameter estimates obtained by the least squares method are 
not unbiased, and as pointed out earlier, the estimates of the 
variances, as developed in this report merely indicate the 
relative accuracy between the different coefficients 

At this point, it seems appropriate to compare the per- 
formance of the two methods. Judging only according to 
efficiency in computer time, the method of quasll Inearization 
is clearly superior to evol ut 1 onary. programml ng for a problem 
with linear differential equations and quadratic cost function. 
In about one minute computer time (CDC 3600) five Iterations 
on a time history with about 150 measured points can be per- 
formed These five iterations yield a set of coefficients 
which minimize the cost function. The coefficients are accur- 
ate to four to five significant figures and as a by-product, 
the estimates of the variances are obtained 

On the other hand, a six minute run on the same computer 
'usinq the evolutionary programming technique yielded a final 
value of the cost function still about twice the size of 
the true mi nimum value Furthermore, no estimates of the error 
bounds are available with this method The evolutionary 
programming technique of minimizing functions with a relative- 
ly large number of unknowns may have advantages in system 
identification problems where a nonlinear model of the plant 
is required or where a nonquadratic error criterion has to 
be used The combination of a general numberlcal integration 
technique (such as for instance Runge-Kutta) and evolutionary 
programming allows quick changes both in the differential 
equations of the model and of the form of the cost function 
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CONCLUSIONS AND RECOMMENDATIONS 


The results of this report show that the method of 
quasi 1 1 near 1 zat i on results in an efficient digital computer 
program which allows determining those values of the stability 
and control derivatives which minimize the integral of the 
weiqhted squared difference between the observed time 
history and the one obtained by solving the differential 
equations using the observed control variables and the para- 
meters to be determined. A few iterations (typically three 
to five) will yield the correct values (the ones which mini- 
mize the cost function) of the unknowns and estimates of their 
variances can be obtained 

The second method which is based on evolutionary program- 
ming cannot compete successfully in the case of linear 
differential equations and a quadratic error function, but it 
may have advantages in nonlinear process identification problems 

The fact that a method is available that solves the 
minimization problem of a given cost function for a given 
form of the systems' differential equations (here linear) 
should not lead to the conclusion that the problem of deter- 
mining stability and control derivatives from flight data is 
solved in its widest engineering sense A number of important 
questions are still open, for instance: 

(1) How does the inclusion or omission of a fit of the 
observed yaw and roll rate derivatives influence the 
coefficients and their variances’ „ 

(2) What criterion should be applied In choosing the 
weighting factors 7 

(3) Is it possible to give variances which are theoreti- 
cally more solidly founded and which distinguish between 
errors due to measurement noise and Inadequacy of the 

- mathematical model 7 

The analysis of the X - 1 5 flight data shown in this report 
seems to indicate that the assumed mathematical model may not 
he quite adequate Especially the fit to the roll rate suggests 
that there are certain terms missing in the roll moment equation, 
these might be unsteady flow derivatives 

The experiments suggest that additional work on the choice 
of the mathematical model with alternative forms of the equations 
of motion (possibly nonlinear equations) be performed Experi- 
ments. where stability and control derivatives obtained from 


52 


one flight test may indicate which mathematical models give 
the most consistent results and therefore are the most 
lealistic ones Computational methods and computer proqrams 
are now available which may help to advance the state of the 
art in the determination and possibly the usage of stability 
and control derivatives 


Decision Science, Inc 
' San Diego, September 25, 1968. 
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